      subroutine ornl_diir2d(    f,incf,ldf,ifx0,n_fx,ify0,n_fy,
     $			    g,incg,ldg,igx0,n_gx,igy0,n_gy,
     $			    h,inch,ldh,ihx0,n_hx,ihy0,n_hy)
      double precision f(ldf,*)
      double precision g(ldg,*)
      double precision h(ldh,*)
      integer incf, ldf, ifx0, n_fx, ify0, n_fy
      integer incg, ldg, igx0, n_gx, igy0, n_gy
      integer inch, ldh, ihx0, n_hx, ihy0, n_hy
c-----------------------------------------------------------------------------
c
c   diir2d  performs a 2D convolution in the time domain :
c	h(i,j) = beta * h(i,j) + alpha * Sum.Sum[ f(k,l) * g(i-k,j-l) ]
c
c-----------------------------------------------------------------------------
c
c   PARAMETERS:
c
c	f:	Vector containing sequence "f"
c	incf:	Increment between two successive lines   of "f"
c	ldf:	Increment between two successive columns of "f"
c	ifx0:	Index of the first element of each column of "f"
c	n_fx:	Number of elements (lines) of each column of "f"
c	ify0:	Index of the first column of "f"
c	n_fy:	Index of the  last column of "f"
c	
c	g:	Vector containing sequence "g"
c	incg:	Increment between two successive lines   of "g"
c	ldg:	Increment between two successive columns of "g"
c	igx0:	Index of the first element of each column of "g"
c	n_gx:	Number of elements (lines) of each column of "g"
c	igy0:	Index of the first column of "g"
c	n_gy:	Index of the  last column of "g"
c	
c	h:	Vector containing sequence "h"
c	inch:	Increment between two successive lines   of "h"
c	ldh:	Increment between two successive columns of "h"
c	ihx0:	Index of the first element of each column of "h"
c	n_hx:	Number of elements (lines) of each column of "h"
c	ihy0:	Index of the first column of "h"
c	n_hy:	Index of the  last column of "h"
c	
c-----------------------------------------------------------------------------
c
c   PLEASE NOTE:
c	Please note that the array pointers must all point to the first 
c	element of the array "(ifx0,ify0)", "(igx0,igy0)" and "(ihx0,ihy0)"
c	 If "f" for example is defined as
c	    dimension f(-25:45,10:21)
c	Then "diir2d" must be called with the following parameters
c	    call diir2d( f(-25,10),(45-(-25)+1),-25,45,10,21 ... )
c
c-----------------------------------------------------------------------------
c
c   This Fortran Subroutine written by
c	Jean-Pierre Panziera
c	Silicon Graphics Inc
c	September 27, 1991
c
c-----------------------------------------------------------------------------

      double precision zero, one
      parameter (zero = 0.0, one = 1.0)
      integer iy0, iy1, ify1, igy1, ihy1
      integer k1, k2, k, j, jf, jg, jh
      integer ix0, ix1, ifx1, igx1, ihx1, if0, ifx, ihx

c-----------------------------------------------------------------------------
	if( (n_hx .le. 0) .or. (n_hy .le. 0) ) return
c-----------------------------------------------------------------------------
c
c   Compute Y Boundaries
c
c-----------------------------------------------------------------------------
	ify1 = ify0 + n_fy - 1
	igy1 = igy0 + n_gy - 1
	ihy1 = ihy0 + n_hy - 1

	iy0 = min(max(ihy0, ify0-igy0), ihy1+1)
	iy1 = max(min(ihy1, ify1-igy0), ihy0-1)

	ifx1 = ifx0 + n_fx - 1
	igx1 = igx0 + n_gx - 1
	ihx1 = ihx0 + n_hx - 1

	ix0 = min( max(ihx0, ifx0-igx0), ihx1+1)
	ix1 = min(ihx1, ifx1-igx0)

c-----------------------------------------------------------------------------
c
c   Zero the Front columns
c
c-----------------------------------------------------------------------------
	do j = ihy0, iy0-1
	    ihx = 1
	    jh = 1 + (j-ihy0)
	    do k = 1, n_hx
		h(ihx,jh) = zero
		ihx = ihx + inch
	    end do
	end do
c-----------------------------------------------------------------------------
c
c   Init the Plain columns
c
c-----------------------------------------------------------------------------
	if0 = 1 + ((ix0+igx0)-ifx0) * incf
	do j = iy0, iy1
	    ihx = 1
	    jh = 1 + (j-ihy0)
	    do k = ihx0, ix0-1
		h(ihx,jh) = zero
		ihx = ihx + inch
	    end do
	    ifx = if0
	    jf = 1 + (j+igy0-ify0)
	    do k = k, ix1
		h(ihx,jh) = f(ifx,jf)
		ihx = ihx + inch
		ifx = ifx + incf
	    end do
	    do k = k, ihx1
		h(ihx,jh) = zero
		ihx = ihx + inch
	    end do
	end do
c-----------------------------------------------------------------------------
c
c   Zero the TAIL columns
c
c-----------------------------------------------------------------------------
	do j = iy1+1, ihy1
	    ihx = 1
	    jh = 1 + (j-ihy0)
	    do k = 1, n_hx
		h(ihx,jh) = zero
		ihx = ihx + inch
	    end do
	end do
c-----------------------------------------------------------------------------
	if( (n_gx .le. 0) .or. (n_gy .le. 0) ) return
c-----------------------------------------------------------------------------
c
c   Call the 1D Convolution
c
c-----------------------------------------------------------------------------
	do j = iy0, ihy1
	    jh = j-ihy0 + 1
	    k1 = max( iy0, j-(n_gy-1))
	    k2 = j-1
	    do k = k1, k2
		jf = k-ihy0 + 1
		jg = j-k + 1
		call dfir1d(	h(1,jf), inch, ihx0, n_hx,
     $				g(1,jg), incg,    0, n_gx,
     $				h(1,jh), inch, ihx0, n_hx,
     $				-one, one)
	    end do
	    call diir1d(	h(1,jh), inch, ihx0, n_hx,
     $				g(1,1 ), incg,    0, n_gx,
     $				h(1,jh), inch, ihx0, n_hx)
	end do
c-----------------------------------------------------------------------------
      return
      end


